Polusi udara merupakan masalah lingkungan yang signifikan di kota-kota padat penduduk, terutama yang menjadi pusat kegiatan ekonomi dan industri, seperti Jakarta, ibu kota Indonesia. Polusi udara di Jakarta sebagian besar berasal dari konsumsi energi, penggunaan bahan bakar fosil, aktivitas rumah tangga, dan proses industri. Di antara polutan udara, partikel materi halus 2.5 (\(PM_{2.5}\)) menjadi ancaman serius bagi kesehatan manusia karena ukurannya yang kecil (≤ 2,5 mikrometer, disingkat µm), yang memungkinkan untuk masuk jauh ke dalam sistem pernapasan. Menyadari bahayanya, Pemerintah Indonesia telah menetapkan ambang batas sehat untuk \(PM_{2.5}\) sebesar 50 µg/m³ (mikrogram per meter kubik). Informasi yang dapat diandalkan tentang konsentrasi \(PM_{2.5}\) sangat penting bagi pembuat kebijakan untuk menyusun strategi mitigasi dan bagi masyarakat untuk melindungi diri dari paparan.
Meskipun penting, pemantauan konsentrasi \(PM_{2.5}\) sering kali terbatas oleh ketersediaan peralatan pemantauan, yang mengakibatkan kekosongan data di berbagai wilayah. Untuk mengatasi hal ini, teknik interpolasi digunakan untuk memperkirakan konsentrasi polutan di lokasi yang tidak terpantau. Metode interpolasi yang umum digunakan meliputi Inverse Distance Weighting (IDW) dan kriging.
IDW memperkirakan nilai di titik yang tidak diketahui hanya berdasarkan jaraknya dari lokasi yang diketahui, dengan memberikan bobot lebih tinggi pada titik yang lebih dekat. Namun, IDW tidak memperhitungkan ketergantungan spasial atau pola variabilitas yang mendasarinya, dan tidak dapat memasukkan variabel prediktor tambahan, yang membatasi kemampuannya untuk memberikan estimasi akurat dalam sistem lingkungan yang kompleks.
Kriging meningkatkan IDW dengan memanfaatkan ketergantungan spasial melalui pemodelan semivariogram, yang menangkap hubungan antara variabilitas dan jarak di antara lokasi yang diamati. Hal ini menghasilkan prediksi yang lebih akurat dibandingkan IDW, terutama di wilayah dengan autokorelasi spasial yang signifikan. Namun, seperti IDW, kriging tidak memungkinkan dimasukkannya variabel prediktor dan tidak secara eksplisit mempertimbangkan ketergantungan temporal.
Namun, polusi udara dipengaruhi oleh kombinasi faktor spasial dan temporal, serta variabel sekunder seperti curah hujan, kelembapan, kecepatan angin, dan konsentrasi \(NO_2\). Untuk mengatasi keterbatasan ini, cokriging spasiotemporal digunakan karena mengintegrasikan ketergantungan spasial dan temporal ke dalam model serta memungkinkan dimasukkannya variabel sekunder sebagai prediktor. Pendekatan ini memungkinkan analisis yang lebih komprehensif terhadap pola distribusi polutan dan meningkatkan akurasi prediksi dengan memasukkan pengaruh faktor lingkungan tambahan bersama konsentrasi \(PM_{2.5}\).
Dalam konteks ini, mengintegrasikan ketergantungan spasial-temporal dan variabel lingkungan tambahan ke dalam model cokriging diharapkan dapat meningkatkan akurasi estimasi konsentrasi \(PM_{2.5}\). Penelitian ini bertujuan untuk mengembangkan model cokriging spasiotemporal untuk memberikan prediksi yang lebih akurat tentang konsentrasi \(PM_{2.5}\) di Jakarta. Dengan memasukkan curah hujan, kelembapan, kecepatan angin, dan konsentrasi \(NO_2\) sebagai variabel sekunder, model yang diusulkan bertujuan untuk menangkap variasi spasial dan temporal secara lebih efektif, sehingga menghasilkan estimasi yang lebih akurat yang mencerminkan dinamika polusi.
library(gstat); library(sf); library(spacetime); library(sp);library(ggplot2)
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(plotly); library(xts); library(RColorBrewer); library(openxlsx)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(raster); library(geojson); library(patchwork); library(ggplotify)
##
## Attaching package: 'raster'
## The following object is masked from 'package:plotly':
##
## select
##
## Attaching package: 'geojson'
## The following object is masked from 'package:graphics':
##
## polygon
##
## Attaching package: 'patchwork'
## The following object is masked from 'package:raster':
##
## area
library(moments); library(tidyr); library(dplyr); library(psych)
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:raster':
##
## extract
##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## ###############################################################################
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:xts':
##
## first, last
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
##
## Attaching package: 'psych'
## The following object is masked from 'package:raster':
##
## distance
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
Pertama, analisis statistik deskriptif dilakukan pada variabel-variabel penelitian.
# INPUT DATA
# PM2.5
pm25 = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx",
sheet = "PM25")
pm25$time = ISOdate(pm25$Tahun, pm25$Bulan, 1, tz = "GMT")
stations = 3:31
pm25_mat = as.matrix(pm25[stations])
# Curah Hujan
presipitasi = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx",
sheet = "Presipitasi")
presipitasi$time = ISOdate(presipitasi$Tahun, presipitasi$Bulan, 1, tz = "GMT")
stations = 3:31
presipitasi_mat = as.matrix(presipitasi[stations])
# Humidity
kelembaban = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx",
sheet = "Kelembaban")
kelembaban$time = ISOdate(kelembaban$Tahun, kelembaban$Bulan, 1, tz = "GMT")
stations = 3:31
kelembaban_mat = as.matrix(kelembaban[stations])
# Data Kecepatan Angin
angin = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx",
sheet = "Kecepatan Angin")
angin$time = ISOdate(angin$Tahun, angin$Bulan, 1, tz = "GMT")
stations = 3:31
angin_mat = as.matrix(angin[stations])
# Data NO2
no2 = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx",
sheet = "NO2")
no2$time = ISOdate(no2$Tahun, no2$Bulan, 1, tz = "GMT")
stations = 3:31
no2_mat = as.matrix(no2[stations])
# MEMBUAT STFDF CLASS OBJECT
tesis.loc = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx",
sheet = "Koordinat")
coordinates(tesis.loc) = ~Longitude+Latitude
proj4string(tesis.loc) = "+proj=longlat +datum=WGS84"
pts = coordinates(tesis.loc)
rownames(pts) = tesis.loc$Sensor
pts = SpatialPoints(pts, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84"))
utm48s = "+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs"
pts.sfc = st_transform(st_as_sfc(pts), utm48s)
pts = as(pts.sfc, "Spatial")
# membuat spatial object dari data masing-masing variabel
pm25.data = stConstruct(pm25_mat, space = list(values = 1:ncol(pm25_mat)),
time = pm25$time, SpatialObj = pts, interval = TRUE)
presipitasi.data = stConstruct(presipitasi_mat, space = list(values = 1:ncol(presipitasi_mat)),
time = presipitasi$time, SpatialObj = pts, interval = TRUE)
kelembaban.data = stConstruct(kelembaban_mat, space = list(values = 1:ncol(kelembaban_mat)),
time = kelembaban$time, SpatialObj = pts, interval = TRUE)
angin.data = stConstruct(angin_mat, space = list(values = 1:ncol(angin_mat)),
time = angin$time, SpatialObj = pts, interval = TRUE)
no2.data = stConstruct(no2_mat, space = list(values = 1:ncol(no2_mat)),
time = no2$time, SpatialObj = pts, interval = TRUE)
# STFDF Object untuk gabungan
combined_data <- cbind(as.data.frame(pm25.data)[,-7],
PM25 = as.data.frame(pm25.data)[,7],
Presipitasi = as.data.frame(presipitasi.data)[,7],
Kelembaban = as.data.frame(kelembaban.data)[,7],
Angin = as.data.frame(angin.data)[,7],
NO2 = as.data.frame(no2.data)[,7])
descriptive_stats <- function(x) {
c(
Min = min(x, na.rm = TRUE),
Max = max(x, na.rm = TRUE),
Mean = mean(x, na.rm = TRUE),
`Std Dev` = sd(x, na.rm = TRUE)
)
}
# Menghitung statistik deskriptif dan membuat hasil dalam format panjang.
results_long <- combined_data[,6:11] %>%
pivot_longer(cols = -timeIndex, names_to = "Variable", values_to = "Value") %>%
group_by(Variable) %>%
summarise(
Minimum = min(Value, na.rm = TRUE),
Maximum = max(Value, na.rm = TRUE),
Mean = mean(Value, na.rm = TRUE),
`Std Dev` = sd(Value, na.rm = TRUE),
.groups = "drop"
)
# Menampilkan hasil
print(results_long)
## # A tibble: 5 × 5
## Variable Minimum Maximum Mean `Std Dev`
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Angin 2.14 6.22 3.79 0.870
## 2 Kelembaban 67.1 89.8 78.4 4.83
## 3 NO2 3.9 30.0 19.4 5.83
## 4 PM25 14.6 68.8 39.6 11.7
## 5 Presipitasi 0 559. 140. 138.
Konsentrasi \(PM_{2.5}\) menunjukkan variasi sepanjang tahun, yang ditunjukkan oleh standar deviasi sebesar 11,69 µg/m³. Nilai ini menunjukkan bahwa konsentrasi \(PM_{2.5}\) menyimpang sebesar 11,69 µg/m³ dari rata-rata. Konsentrasi \(PM_{2.5}\) terendah terjadi pada bulan Februari, dengan nilai 14,56 µg/m³, sedangkan konsentrasi maksimum tercatat pada bulan Agustus, mencapai 68,83 µg/m³. Konsentrasi \(NO_{2}\) terendah adalah 3,9 µg/m³, yang terjadi pada bulan Desember, sedangkan konsentrasi tertinggi, juga tercatat pada bulan Desember tetapi di lokasi yang berbeda, adalah 29,97 µg/m³. Curah hujan menunjukkan variabilitas yang lebih tinggi dibandingkan dengan variabel lainnya, sebagaimana tercermin dalam standar deviasi sebesar 137,78 mm. Curah hujan minimum di Jakarta selama periode 2023 adalah 0 mm, yang menunjukkan tidak ada hujan sama sekali dalam satu bulan penuh, sedangkan curah hujan maksimum mencapai 558,83 mm. Variabel kelembapan udara dan kecepatan angin menunjukkan nilai yang relatif homogen di seluruh lokasi dan periode waktu, dengan standar deviasi masing-masing sebesar 4,82% dan 0,87 m/s.
Terdapat korelasi negatif antara curah hujan, kelembapan udara, dan kecepatan angin terhadap konsentrasi \(PM_{2.5}\). Ketika kelembapan udara, kecepatan angin, dan curah hujan meningkat, konsentrasi \(PM_{2.5}\) cenderung menurun. Sebaliknya, \(NO_{2}\) memiliki korelasi positif dengan \(PM_{2.5}\). Ini berarti bahwa konsentrasi \(PM_{2.5}\) juga akan meningkat seiring dengan meningkatnya \(NO_{2}\).
corr.test(combined_data[,c(7,11,8:10)], method = 'pearson', adjust = 'none')
## Call:corr.test(x = combined_data[, c(7, 11, 8:10)], method = "pearson",
## adjust = "none")
## Correlation matrix
## PM25 NO2 Presipitasi Kelembaban Angin
## PM25 1.00 0.31 -0.58 -0.13 -0.75
## NO2 0.31 1.00 -0.34 -0.43 -0.19
## Presipitasi -0.58 -0.34 1.00 0.73 0.54
## Kelembaban -0.13 -0.43 0.73 1.00 0.04
## Angin -0.75 -0.19 0.54 0.04 1.00
## Sample Size
## [1] 348
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## PM25 NO2 Presipitasi Kelembaban Angin
## PM25 0.00 0 0 0.02 0.00
## NO2 0.00 0 0 0.00 0.00
## Presipitasi 0.00 0 0 0.00 0.00
## Kelembaban 0.02 0 0 0.00 0.47
## Angin 0.00 0 0 0.47 0.00
##
## To see confidence intervals of the correlations, print with the short=FALSE option
Dari total sampel sebanyak 29 titik, data akan dibagi menjadi 80% data pelatihan (23 titik) dan 20% data pengujian (6 titik).
# Baca data lokasi stasiun dari file Excel
sensor_loc <- read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx",
sheet = 'Koordinat')
# RANDOM SAMPLE
sample_index = unique(combined_data$sp.ID)
#set.seed(123) # untuk konsistensi
#test_indices <- sample(sample_index, size = round(0.2 * length(sample_index),0))
#saveRDS(test_indices,
# "D:/TESIS/DATA/Data Fix/indeks wilayah test.RDS")
test_samples = readRDS("D:/TESIS/DATA/Data Fix/indeks wilayah test.RDS")
train_samples = as.numeric(setdiff(sample_index, test_samples))
# Convert the data frame to a spatial points data frame
stasiun_test = st_as_sf(sensor_loc[sensor_loc$Sensor %in% c("Cibubur", "Cipinang Besar",
"Duri Utara", "Marunda",
"Rempoa Permai", "Semanan"),],
coords = c("Longitude", "Latitude"), crs = 4326)
stasiun_train = st_as_sf(sensor_loc[!(sensor_loc$Sensor %in% c("Cibubur", "Cipinang Besar",
"Duri Utara", "Marunda",
"Rempoa Permai", "Semanan")),],
coords = c("Longitude", "Latitude"), crs = 4326)
# Load the Jakarta boundary shapefile (adjust the path to your shapefile)
province_boundary = st_read("D:/TESIS/DATA/SHP/gadm41_IDN_2.shp")
## Reading layer `gadm41_IDN_2' from data source `D:\TESIS\DATA\SHP\gadm41_IDN_2.shp' using driver `ESRI Shapefile'
## Simple feature collection with 502 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.00971 ymin: -11.00761 xmax: 141.0194 ymax: 6.076941
## Geodetic CRS: WGS 84
# Inspect the structure of the province_boundary
str(province_boundary)
## Classes 'sf' and 'data.frame': 502 obs. of 14 variables:
## $ GID_2 : chr "IDN.1.2_1" "IDN.1.1_1" "IDN.1.3_1" "IDN.1.4_1" ...
## $ GID_0 : chr "IDN" "IDN" "IDN" "IDN" ...
## $ COUNTRY : chr "Indonesia" "Indonesia" "Indonesia" "Indonesia" ...
## $ GID_1 : chr "IDN.1_1" "IDN.1_1" "IDN.1_1" "IDN.1_1" ...
## $ NAME_1 : chr "Aceh" "Aceh" "Aceh" "Aceh" ...
## $ NL_NAME_1: chr "NA" "NA" "NA" "NA" ...
## $ NAME_2 : chr "Aceh Barat" "Aceh Barat Daya" "Aceh Besar" "Aceh Jaya" ...
## $ VARNAME_2: chr "NA" "NA" "NA" "NA" ...
## $ NL_NAME_2: chr "NA" "NA" "NA" "NA" ...
## $ TYPE_2 : chr "Kabupaten" "Kabupaten" "Kabupaten" "Kabupaten" ...
## $ ENGTYPE_2: chr "Regency" "Regency" "Regency" "Regency" ...
## $ CC_2 : chr "1107" "1112" "1108" "1116" ...
## $ HASC_2 : chr "ID.AC.AB" "ID.AC.AD" "ID.AC.AR" "ID.AC.AJ" ...
## $ geometry :sfc_MULTIPOLYGON of length 502; first list element: List of 1
## ..$ :List of 1
## .. ..$ : num [1:929, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
## ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:13] "GID_2" "GID_0" "COUNTRY" "GID_1" ...
jkt_boundary <- province_boundary %>%
filter(NAME_1 == "Jakarta Raya")
#jkt_boundary[-6,"NAME_EN"] = c('West Jakarta', 'Central Jakarta', 'South Jakarta',
# 'East Jakarta', 'North Jakarta')
jkt_boundary[-6,"NAME_EN"] = c('Jakarta Barat', 'Jakarta Pusat', 'Jakarta Selatan',
'Jakarta Timur', 'Jakarta Utara')
# Extract the boundary lines for visualization
boundary_lines <- st_cast(jkt_boundary, "MULTILINESTRING")
# Plot the points and the filtered Jakarta boundary with legend title as "SPKU"
ggplot() +
geom_sf(data = jkt_boundary[-6,], fill = "lightyellow", color = "brown") +
geom_sf_text(data = jkt_boundary[-6,], aes(label = NAME_EN), size = 4,
nudge_y = c(0.01, 0.01, 0.02, 0.02, 0.01),
nudge_x = c(0, 0.02, 0, 0.02, 0), color = "black") +
geom_sf(data = stasiun_train, aes(color = "Latih"), size = 2) + # Tambahkan label legenda
geom_sf(data = stasiun_test, aes(color = "Uji"), size = 2) + # Tambahkan label legenda
scale_color_manual(name = "Data", values = c("Latih" = "blue", "Uji" = "orangered")) + # Skala warna
theme_minimal() +
labs(x = "Longitude",
y = "Latitude") +
theme(panel.grid.major = element_blank(), # Ubah grid besar jadi putih
panel.grid.minor = element_blank(), # Ubah grid kecil jadi putih
legend.title = element_text(size = 12), # Ukuran judul legenda
legend.text = element_text(size = 10)) # Ukuran teks legenda
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
Gambar tersebut menunjukkan pembagian data yang digunakan sebagai data latih dan data uji. Observasi yang ditunjukkan oleh titik berwarna biru menunjukkan data latih, sedangkan titik berwarna merah menunjukkan data uji. Berikut ini adalah proses membagi data dan mengubah data ke dalam format STFDF
# Konversi kolom waktu ke format POSIXct
Sys.setenv(TZ = "GMT")
time = as.POSIXct(unique(combined_data$time), tz = "GMT")
# RANDOM SAMPLE
sample_index = unique(combined_data$sp.ID)
#set.seed(123) # untuk konsistensi
#test_indices <- sample(sample_index, size = round(0.2 * length(sample_index),0))
#saveRDS(test_indices,
# "D:/TESIS/DATA/Data Fix/indeks wilayah test.RDS")
test_samples = readRDS("D:/TESIS/DATA/Data Fix/indeks wilayah test.RDS")
# Membuat objek STFDF DATA TRAIN
train_samples = as.numeric(setdiff(sample_index, test_samples))
train_data = combined_data[combined_data$sp.ID %in% train_samples, ]
coords_train = unique(train_data[, c("coords.x1","coords.x2")])
pts_train = coordinates(coords_train)
pts_train = SpatialPoints(pts_train,
CRS("+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs"))
utm48s = "+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs"
pts_train.sfc = st_transform(st_as_sfc(pts_train), utm48s)
pts_train = as(pts_train.sfc, "Spatial")
stfdf.train = STFDF(sp = pts_train, time = time, data = train_data[,7:11])
stfdf.train
## An object of class "STFDF"
## Slot "data":
## PM25 Presipitasi Kelembaban Angin NO2
## 1 25.26000 169.67 80.50000 4.87 9.37000
## 2 29.50000 185.37 82.27000 4.64 7.97000
## 3 24.23000 165.96 80.07000 4.92 9.71000
## 5 31.47000 192.90 83.10000 4.53 7.31000
## 6 23.65000 163.89 79.83000 4.95 9.90000
## 8 28.92000 183.18 82.03000 4.67 8.16000
## 10 21.32000 155.70 78.85000 5.08 10.67000
## 11 21.18000 155.22 78.79000 5.09 10.71000
## 12 23.93667 164.91 79.95000 4.94 9.80000
## 13 26.77000 175.18 81.13000 4.78 8.87000
## 14 20.47000 152.77 78.50000 5.13 21.60000
## 15 21.39000 155.95 78.88000 5.08 10.64000
## 16 20.42000 152.60 78.48000 5.13 10.96000
## 17 28.30000 180.86 81.77000 4.70 8.36000
## 19 22.67000 160.42 79.42000 5.01 10.22000
## 20 17.01000 141.11 77.05000 5.31 12.09000
## 21 29.31000 184.65 82.19000 4.65 8.03000
## 22 22.77000 160.77 79.46000 5.00 10.19000
## 23 22.07000 158.32 79.17000 5.04 10.42000
## 26 18.78000 147.02 77.79000 5.22 11.51000
## 27 21.62000 169.57 85.69774 5.40 10.57000
## 28 22.68000 160.46 79.42000 5.01 10.22000
## 29 22.80000 160.88 79.47000 5.00 10.18000
## 30 20.87000 502.43 84.09000 5.88 14.36000
## 31 23.76000 520.60 85.30000 5.72 13.41000
## 32 21.88000 508.74 84.52000 5.82 14.03000
## 34 23.49000 518.89 85.19000 5.73 13.50000
## 35 20.35000 499.20 83.88000 5.90 14.54000
## 37 24.85000 527.54 85.76000 5.66 13.05000
## 39 19.31333 492.78 83.44000 5.96 14.88000
## 40 21.85000 508.56 84.50000 5.82 14.04000
## 41 21.23333 504.70 84.25000 5.86 14.24000
## 42 23.01000 515.85 84.99000 5.76 13.66000
## 43 17.78000 483.36 82.80000 6.04 18.64286
## 44 20.34000 499.13 83.87000 5.90 14.54000
## 45 17.92000 484.22 82.86000 6.04 15.34000
## 46 23.57000 519.40 85.22000 5.73 13.47000
## 48 20.51000 500.19 83.94000 5.89 14.48000
## 49 14.56000 463.88 81.46000 6.22 16.45000
## 50 26.45000 537.80 86.43000 5.57 12.52000
## 51 18.59000 488.32 83.14000 6.00 15.12000
## 52 19.28000 492.57 83.43000 5.96 14.89000
## 55 20.15000 497.95 83.79000 5.91 14.60000
## 56 18.82000 442.02 86.65000 5.60 15.04000
## 57 20.23000 498.45 83.83000 5.91 14.58000
## 58 21.18000 504.36 84.22000 5.86 14.26000
## 59 39.89000 285.12 86.16000 3.54 20.24000
## 60 38.60000 279.10 85.62000 3.61 20.66000
## 61 31.84000 248.58 82.79000 3.97 22.89000
## 63 42.27000 296.40 87.16000 3.41 19.45000
## 64 32.73000 252.50 83.17000 3.92 22.60000
## 66 40.03000 285.78 86.22000 3.53 20.19000
## 68 30.09667 241.00 82.07000 4.07 23.47000
## 69 32.92000 253.34 83.25000 3.91 22.54000
## 70 29.21000 237.19 81.69000 4.12 23.76000
## 71 39.38000 282.73 85.95000 3.56 20.40000
## 72 31.91000 248.89 82.82000 3.97 24.80645
## 73 31.48000 247.01 82.64000 3.99 23.01000
## 74 28.74000 235.18 81.50000 4.14 23.92000
## 75 37.90000 275.85 85.33000 3.64 20.89000
## 77 29.69000 239.25 81.90000 4.09 23.60000
## 78 23.83000 214.72 79.45000 4.41 25.54000
## 79 37.80000 275.39 85.29000 3.65 20.93000
## 80 31.27000 246.09 82.56000 4.00 23.08000
## 81 26.89000 227.36 80.72000 4.24 24.53000
## 84 33.18000 254.49 83.35000 3.90 22.45000
## 85 27.82000 231.10 84.35871 3.40 24.22000
## 86 28.18000 232.80 81.26000 4.17 24.10000
## 87 28.65000 234.80 81.46000 4.15 23.95000
## 88 32.54000 182.31 80.93000 3.47 15.78000
## 89 34.40000 189.36 81.71000 3.37 15.17000
## 90 25.49000 156.81 77.99000 3.85 18.11000
## 92 37.67000 202.08 83.08000 3.19 14.09000
## 93 28.02000 165.74 79.04000 3.72 17.28000
## 95 34.63000 190.24 81.81000 3.36 15.09000
## 97 24.70000 154.07 77.66000 3.90 18.37000
## 98 26.77000 161.30 78.52000 3.78 17.69000
## 99 25.49000 156.81 77.99000 3.85 18.11000
## 100 33.20000 184.80 81.21000 3.44 15.56000
## 101 28.21000 166.42 79.12000 3.71 18.86667
## 102 23.85000 151.15 77.30000 3.94 18.65000
## 103 25.16000 155.66 77.85000 3.87 18.22000
## 104 32.02000 180.37 80.72000 3.50 15.95000
## 106 25.87000 158.14 78.15000 3.83 17.99000
## 107 22.79000 147.55 76.86000 4.00 19.00000
## 108 33.57000 186.20 81.37000 3.42 15.44000
## 109 28.30000 166.74 79.16000 3.70 17.18000
## 110 23.42000 149.69 77.12000 3.97 18.79000
## 113 29.89000 172.50 79.83000 3.61 16.66000
## 114 26.12000 114.66 83.01833 4.50 17.90000
## 115 25.98000 158.52 78.19000 3.83 17.95000
## 116 26.13000 159.05 78.25000 3.82 17.90000
## 117 55.43000 164.05 84.00000 2.72 20.84000
## 118 58.48000 175.09 85.27000 2.56 19.84000
## 119 36.25000 102.87 75.98000 3.76 27.18000
## 121 60.11000 181.14 85.95000 2.47 19.30000
## 122 44.82000 128.45 79.56000 3.30 24.35000
## 124 53.24000 156.35 83.08000 2.84 21.57000
## 126 49.44000 143.42 81.49000 3.05 22.82000
## 127 46.85000 134.92 80.41000 3.19 23.68000
## 128 39.72000 112.88 77.43000 3.57 26.03000
## 129 53.34000 156.69 83.12000 2.83 21.53000
## 130 36.65000 104.00 76.14000 3.74 22.56667
## 131 42.10000 120.02 78.42000 3.44 25.25000
## 132 41.65000 118.66 78.23000 3.47 25.40000
## 133 50.75000 147.81 82.04000 2.98 22.39000
## 135 40.28000 114.54 77.66000 3.54 25.85000
## 136 37.94000 107.69 76.68000 3.67 26.62000
## 137 49.91000 144.99 81.69000 3.02 22.67000
## 138 46.94000 135.22 80.45000 3.18 23.65000
## 139 34.64000 98.38 75.30000 3.85 27.71000
## 142 46.32000 133.22 80.19000 3.22 23.85000
## 143 42.95000 172.48 81.13839 3.80 24.97000
## 144 32.01000 91.26 74.20000 3.99 28.58000
## 145 40.56000 115.38 77.78000 3.53 25.76000
## 146 53.91000 161.11 84.71000 2.56 17.82000
## 147 49.31000 145.29 82.79000 2.81 19.34000
## 148 38.25000 110.59 78.17000 3.41 22.99000
## 150 53.99000 161.39 84.75000 2.56 17.79000
## 151 41.67000 120.81 79.60000 3.22 21.86000
## 153 51.57000 152.96 83.74000 2.69 18.59000
## 155 42.84000 124.42 80.09000 3.16 21.48000
## 156 42.24000 122.56 79.83000 3.19 21.68000
## 157 40.42000 117.02 79.07000 3.29 22.28000
## 158 49.76000 146.80 82.98000 2.78 19.19000
## 159 43.93000 127.82 80.54000 3.10 22.13333
## 160 41.61000 120.63 79.57000 3.23 21.88000
## 161 40.73000 117.96 79.20000 3.27 22.17000
## 162 48.47000 142.48 82.44000 2.85 19.62000
## 164 39.09000 113.06 78.52000 3.36 22.72000
## 165 35.60000 102.98 77.06000 3.55 23.87000
## 166 50.86000 150.53 83.44000 2.73 18.83000
## 167 45.93000 134.18 81.38000 2.99 20.46000
## 168 34.65000 100.32 76.66000 3.60 24.18000
## 171 42.48000 123.30 79.94000 3.18 21.60000
## 172 43.32000 107.15 82.91300 3.70 21.32000
## 173 36.72000 106.16 77.53000 3.49 23.50000
## 174 40.03000 115.85 78.91000 3.31 22.41000
## 175 56.96000 27.28 80.30000 3.02 20.62000
## 176 51.91000 20.44 78.19000 3.29 22.28000
## 177 38.28000 6.90 72.49000 4.03 26.79000
## 179 67.43000 44.59 84.68000 2.45 17.16000
## 180 47.15000 14.89 76.20000 3.55 23.86000
## 182 55.72000 25.51 79.79000 3.09 21.03000
## 184 41.39667 9.36 73.80000 3.86 25.76000
## 185 48.07000 15.90 76.59000 3.50 23.55000
## 186 46.51000 14.21 75.93000 3.58 24.07000
## 187 50.47000 18.67 77.59000 3.37 22.76000
## 188 51.38000 19.78 77.97000 3.32 24.83871
## 189 43.55000 11.28 74.70000 3.74 25.05000
## 190 44.01000 11.72 74.89000 3.72 24.89000
## 191 54.17000 23.38 79.14000 3.17 21.54000
## 193 43.01000 10.78 74.47000 3.77 25.22000
## 194 38.96000 7.40 72.78000 3.99 26.56000
## 195 55.94000 25.82 79.88000 3.07 20.95000
## 196 46.85000 14.57 76.08000 3.57 23.96000
## 197 37.95000 6.66 72.35000 4.05 26.90000
## 200 48.03000 15.85 76.57000 3.50 23.57000
## 201 47.11000 38.60 79.39516 3.50 23.87000
## 202 47.64000 15.42 76.41000 3.52 23.70000
## 203 44.08000 11.78 74.92000 3.72 24.87000
## 204 56.35000 13.34 78.64000 3.31 19.90000
## 205 51.58000 8.94 76.65000 3.57 21.48000
## 206 40.50000 2.10 72.02000 4.17 25.14000
## 208 67.73000 27.40 83.40000 2.69 16.15000
## 209 50.43000 8.01 76.17000 3.63 21.86000
## 211 54.67000 11.69 77.94000 3.40 20.46000
## 213 41.52000 2.53 72.44000 4.11 24.80000
## 214 48.20000 6.35 75.24000 3.75 22.60000
## 215 48.61000 6.64 75.41000 3.73 22.46000
## 216 49.93000 7.62 75.96000 3.66 22.03000
## 217 49.22000 7.08 75.66000 3.70 19.38710
## 218 45.59000 4.65 74.14000 3.89 23.46000
## 219 42.18000 2.83 72.72000 4.08 24.59000
## 220 56.43000 13.42 78.68000 3.31 19.88000
## 222 42.00000 2.75 72.64000 4.09 24.65000
## 223 40.50000 2.10 72.02000 4.17 25.14000
## 224 56.48000 13.47 78.70000 3.30 19.86000
## 225 47.53000 5.89 74.96000 3.79 22.82000
## 226 37.54000 1.08 70.78000 4.33 26.12000
## 229 50.50000 8.06 76.20000 3.63 21.84000
## 230 48.40000 17.49 78.45774 4.80 22.53000
## 231 43.46000 3.46 73.25000 4.01 24.16000
## 232 44.52000 4.03 73.70000 3.95 23.81000
## 233 51.14000 3.68 73.65000 3.51 21.21000
## 234 44.86000 1.09 71.03000 3.85 23.28000
## 235 39.24000 0.07 68.68000 4.16 25.14000
## 237 60.00000 9.91 77.36000 3.03 18.28000
## 238 42.91000 0.60 70.21000 3.96 23.92000
## 240 48.10000 2.23 72.38000 3.68 22.21000
## 242 41.36000 0.31 69.56000 4.04 24.44000
## 243 44.19000 0.90 70.75000 3.89 23.50000
## 244 42.80000 0.57 70.17000 3.96 23.96000
## 245 42.07000 0.43 69.86000 4.00 24.20000
## 246 43.47000 0.72 70.45000 3.93 20.44828
## 247 42.40000 0.49 70.00000 3.99 24.09000
## 248 38.36000 0.02 68.31000 4.20 25.43000
## 249 46.79000 1.72 71.84000 3.75 22.64000
## 251 35.40000 0.07 67.07000 4.36 26.41000
## 252 37.22000 0.00 67.83000 4.27 25.80000
## 253 53.40000 4.98 74.60000 3.39 20.46000
## 254 41.96000 0.41 69.82000 4.01 24.24000
## 255 37.02000 0.00 67.75000 4.28 25.87000
## 258 44.12000 0.89 70.72000 3.89 23.52000
## 259 46.53000 2.17 73.03300 4.30 22.73000
## 260 44.43000 0.97 70.85000 3.88 23.42000
## 261 42.35000 0.48 69.98000 3.99 24.11000
## 262 52.99000 9.57 73.18000 3.69 21.90000
## 263 58.55000 14.95 75.50000 3.38 20.07000
## 264 41.38000 2.19 68.32000 4.32 25.74000
## 266 62.27000 19.22 77.06000 3.18 18.84000
## 267 48.92000 6.39 71.48000 3.91 23.25000
## 269 52.56000 9.20 73.00000 3.71 22.04000
## 271 46.94000 5.07 70.65000 4.01 23.90000
## 272 50.34000 7.43 72.07000 3.83 22.78000
## 273 44.33000 3.57 69.56000 4.16 24.76000
## 274 48.69000 6.23 71.38000 3.92 23.32000
## 275 48.90000 6.38 71.47000 3.91 21.38710
## 276 45.97000 4.49 70.24000 4.07 24.22000
## 277 40.61000 1.88 68.00000 4.36 25.99000
## 278 53.80000 10.28 73.52000 3.64 21.64000
## 280 41.85000 2.39 68.52000 4.29 25.58000
## 281 40.48000 1.84 67.95000 4.36 26.04000
## 282 59.24000 15.70 75.79000 3.35 19.84000
## 283 46.57000 4.85 70.49000 4.03 24.02000
## 284 40.12000 1.70 67.80000 4.38 26.15000
## 287 47.90000 5.69 71.05000 3.96 23.58000
## 288 51.19000 34.50 71.84032 4.30 22.50000
## 289 48.43000 6.05 71.27000 3.93 23.41000
## 290 47.35000 5.33 70.82000 3.99 23.77000
## 291 42.65000 188.29 78.17000 3.08 23.90000
## 292 53.91000 233.69 82.88000 2.47 20.18000
## 293 36.90000 167.00 75.77000 3.39 25.80000
## 295 47.24000 206.21 80.09000 2.83 22.39000
## 296 40.11000 178.73 77.11000 3.21 24.74000
## 298 45.07000 197.64 79.18000 2.94 23.10000
## 300 40.48000 180.11 77.26000 3.19 24.62000
## 301 43.06000 189.86 78.34000 3.05 23.77000
## 302 36.82000 166.71 75.73000 3.39 25.83000
## 303 38.72000 173.60 76.53000 3.29 25.20000
## 304 42.66000 188.33 78.17000 3.08 29.62069
## 305 39.48000 176.40 76.84000 3.25 24.95000
## 306 30.90000 146.14 73.26000 3.71 27.79000
## 307 44.30000 194.64 78.86000 2.99 23.36000
## 309 35.58000 162.29 75.21000 3.46 26.24000
## 310 34.72000 159.26 74.85000 3.51 26.52000
## 311 51.31000 222.78 81.79000 2.61 21.04000
## 312 38.91000 174.30 76.61000 3.28 25.14000
## 313 35.04000 160.38 74.99000 3.49 26.42000
## 316 40.16000 178.92 77.13000 3.21 24.73000
## 317 41.58000 166.61 76.67733 3.40 24.26000
## 318 37.19000 168.04 75.89000 3.37 25.71000
## 319 41.50000 183.93 77.69000 3.14 24.28000
## 320 40.98000 89.25 77.99000 2.85 7.37000
## 321 45.70000 102.08 79.96000 2.59 5.81000
## 322 32.53000 68.44 74.46000 3.31 10.16000
## 324 39.30000 84.90 77.29000 2.94 7.92000
## 325 36.81000 78.64 76.25000 3.07 8.75000
## 327 42.16000 92.38 78.48000 2.78 6.98000
## 329 35.83000 76.24 75.84000 3.13 9.07000
## 330 38.32000 82.40 76.88000 2.99 8.25000
## 331 33.19000 69.97 74.73000 3.27 9.94000
## 332 42.57000 93.48 78.65000 2.76 6.84000
## 333 38.29000 82.33 76.87000 2.99 29.96774
## 334 33.86000 71.53 75.01000 3.23 9.72000
## 335 27.71000 57.81 72.44000 3.57 11.75000
## 336 41.84000 91.53 78.35000 2.80 7.08000
## 338 31.42000 65.91 73.99000 3.37 10.53000
## 339 30.78000 64.48 73.72000 3.40 10.74000
## 340 46.12000 103.26 80.14000 2.57 5.67000
## 341 35.08000 74.43 75.52000 3.17 9.32000
## 342 30.89000 64.72 73.77000 3.40 10.70000
## 345 35.72000 75.97 75.79000 3.13 9.11000
## 346 36.48000 128.27 77.53710 3.40 8.86000
## 347 30.72000 64.34 73.70000 3.40 10.76000
## 348 37.92000 81.40 76.71000 3.01 8.38000
##
## Slot "sp":
## SpatialPoints:
## coords.x1 coords.x2
## [1,] 701078.1 9317300
## [2,] 702363.8 9315200
## [3,] 703520.7 9314565
## [4,] 702072.5 9315722
## [5,] 713562.7 9314648
## [6,] 708202.7 9297513
## [7,] 708365.4 9311897
## [8,] 705238.6 9305209
## [9,] 711732.0 9317081
## [10,] 711496.6 9321276
## [11,] 716938.0 9325901
## [12,] 712042.2 9316295
## [13,] 697455.1 9322826
## [14,] 706048.1 9320698
## [15,] 700097.1 9312272
## [16,] 702312.6 9304236
## [17,] 696675.0 9302755
## [18,] 694931.3 9305603
## [19,] 703922.7 9302476
## [20,] 693092.9 9317830
## [21,] 697690.2 9314998
## [22,] 688138.4 9317314
## [23,] 698100.5 9317878
## Coordinate Reference System (CRS) arguments: +proj=utm +zone=48 +south
## +datum=WGS84 +units=m +no_defs
##
## Slot "time":
## timeIndex
## 2023-01-01 12:00:00 1
## 2023-02-01 12:00:00 2
## 2023-03-01 12:00:00 3
## 2023-04-01 12:00:00 4
## 2023-05-01 12:00:00 5
## 2023-06-01 12:00:00 6
## 2023-07-01 12:00:00 7
## 2023-08-01 12:00:00 8
## 2023-09-01 12:00:00 9
## 2023-10-01 12:00:00 10
## 2023-11-01 12:00:00 11
## 2023-12-01 12:00:00 12
##
## Slot "endTime":
## [1] "2023-02-01 12:00:00 GMT" "2023-03-01 12:00:00 GMT"
## [3] "2023-04-01 12:00:00 GMT" "2023-05-01 12:00:00 GMT"
## [5] "2023-06-01 12:00:00 GMT" "2023-07-01 12:00:00 GMT"
## [7] "2023-08-01 12:00:00 GMT" "2023-09-01 12:00:00 GMT"
## [9] "2023-10-01 12:00:00 GMT" "2023-11-01 12:00:00 GMT"
## [11] "2023-12-01 12:00:00 GMT" "2023-12-31 12:00:00 GMT"
# STFDF DATA TEST
test_data = combined_data[combined_data$sp.ID %in% test_samples, ]
coords_test = unique(test_data[, c("coords.x1","coords.x2")])
pts_test = coordinates(coords_test)
pts_test = SpatialPoints(pts_test,
CRS("+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs"))
utm48s = "+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs"
pts_test.sfc = st_transform(st_as_sfc(pts_test), utm48s)
pts_test = as(pts_test.sfc, "Spatial")
stfdf.test = STFDF(sp = pts_test, time = time, data = test_data[,7:11])
stfdf.test
## An object of class "STFDF"
## Slot "data":
## PM25 Presipitasi Kelembaban Angin NO2
## 4 38.08 219.27 85.86 4.17 5.13
## 7 18.34 145.54 77.61 5.24 11.65
## 9 27.43 177.62 81.41 4.75 8.65
## 18 22.42 159.54 79.31 5.02 10.30
## 24 36.73 213.74 85.30 4.25 5.58
## 25 30.50 189.17 82.69 4.58 7.63
## 33 28.26 549.53 87.18 5.48 11.92
## 36 15.72 470.85 81.94 6.15 16.07
## 38 24.23 523.59 85.50 5.69 13.25
## 47 18.64 488.63 83.16 6.00 15.10
## 53 29.68 558.83 87.78 5.40 11.45
## 54 29.52 557.78 87.71 5.41 11.51
## 62 46.60 317.48 88.97 3.17 18.02
## 65 26.97 227.70 80.76 4.24 24.50
## 67 29.17 237.02 81.68 4.12 23.78
## 76 29.54 238.61 81.83 4.10 23.65
## 82 46.47 316.84 88.91 3.18 18.06
## 83 48.57 327.32 89.79 3.07 17.37
## 91 41.02 215.53 84.48 3.01 12.98
## 94 24.76 154.28 77.68 3.89 18.35
## 96 26.41 160.03 78.37 3.80 17.81
## 105 30.35 174.18 80.02 3.59 16.51
## 111 35.21 192.47 82.05 3.33 14.90
## 112 34.84 191.05 81.90 3.35 15.02
## 120 64.90 199.50 87.96 2.21 17.72
## 123 40.47 115.11 77.74 3.53 25.79
## 125 37.37 106.05 76.44 3.70 26.81
## 134 50.78 147.91 82.05 2.97 22.38
## 140 58.24 174.21 85.17 2.57 19.92
## 141 41.48 118.14 78.16 3.48 25.45
## 149 61.65 189.58 87.95 2.14 15.26
## 152 39.84 115.29 78.83 3.32 22.47
## 154 36.99 106.94 77.64 3.48 23.41
## 163 49.96 147.47 83.06 2.77 19.13
## 169 53.64 160.16 84.60 2.57 17.91
## 170 58.95 179.38 86.82 2.29 16.16
## 178 64.58 39.46 83.49 2.61 18.10
## 181 42.36 10.20 74.20 3.81 25.44
## 183 41.01 9.03 73.63 3.88 25.89
## 192 53.05 21.90 78.67 3.23 21.91
## 198 56.32 26.36 80.04 3.05 20.83
## 199 64.97 40.15 83.65 2.58 17.97
## 207 68.83 29.02 83.86 2.63 15.78
## 210 42.95 3.20 73.04 4.04 24.33
## 212 40.98 2.30 72.22 4.14 24.98
## 221 52.30 9.54 76.95 3.53 21.24
## 227 58.65 15.78 79.61 3.19 19.15
## 228 67.18 26.60 83.17 2.72 16.33
## 236 62.13 11.87 78.25 2.92 17.58
## 239 41.46 0.33 69.61 4.04 24.40
## 241 39.94 0.13 68.97 4.12 24.91
## 250 48.13 2.25 72.40 3.67 22.20
## 256 53.15 4.83 74.49 3.40 20.54
## 257 62.53 12.25 78.42 2.89 17.44
## 265 65.71 23.64 78.50 3.00 17.70
## 268 41.12 2.08 68.21 4.33 25.82
## 270 44.18 3.49 69.49 4.16 24.81
## 279 54.72 11.12 73.90 3.59 21.33
## 285 61.22 17.96 76.62 3.24 19.18
## 286 68.15 27.05 79.52 2.86 16.89
## 294 54.38 235.70 83.07 2.44 20.03
## 297 33.43 154.77 74.31 3.58 26.95
## 299 38.03 171.08 76.24 3.33 25.43
## 308 44.16 194.10 78.80 2.99 23.40
## 314 53.41 231.57 82.67 2.49 20.35
## 315 54.37 235.65 83.07 2.44 20.03
## 323 51.49 118.99 82.38 2.28 3.90
## 326 30.51 63.87 73.61 3.42 10.83
## 328 34.24 72.43 75.17 3.21 9.60
## 337 40.02 86.75 77.59 2.90 7.69
## 343 45.79 102.33 80.00 2.59 5.78
## 344 47.62 107.54 80.77 2.49 5.18
##
## Slot "sp":
## SpatialPoints:
## coords.x1 coords.x2
## [1,] 702838.7 9315468
## [2,] 706282.9 9308147
## [3,] 712173.4 9301736
## [4,] 699318.8 9306375
## [5,] 690991.0 9322319
## [6,] 699824.2 9319688
## Coordinate Reference System (CRS) arguments: +proj=utm +zone=48 +south
## +datum=WGS84 +units=m +no_defs
##
## Slot "time":
## timeIndex
## 2023-01-01 12:00:00 1
## 2023-02-01 12:00:00 2
## 2023-03-01 12:00:00 3
## 2023-04-01 12:00:00 4
## 2023-05-01 12:00:00 5
## 2023-06-01 12:00:00 6
## 2023-07-01 12:00:00 7
## 2023-08-01 12:00:00 8
## 2023-09-01 12:00:00 9
## 2023-10-01 12:00:00 10
## 2023-11-01 12:00:00 11
## 2023-12-01 12:00:00 12
##
## Slot "endTime":
## [1] "2023-02-01 12:00:00 GMT" "2023-03-01 12:00:00 GMT"
## [3] "2023-04-01 12:00:00 GMT" "2023-05-01 12:00:00 GMT"
## [5] "2023-06-01 12:00:00 GMT" "2023-07-01 12:00:00 GMT"
## [7] "2023-08-01 12:00:00 GMT" "2023-09-01 12:00:00 GMT"
## [9] "2023-10-01 12:00:00 GMT" "2023-11-01 12:00:00 GMT"
## [11] "2023-12-01 12:00:00 GMT" "2023-12-31 12:00:00 GMT"
# membuat object gstat untuk variogram
g <- gstat(NULL, "PM25", PM25 ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
gstat("Presipitasi", Presipitasi ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
gstat("Kelembaban", Kelembaban ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
gstat("Angin", Angin ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
gstat("NO2", NO2 ~ 1, data = stfdf.train, fill.cross = TRUE)
# =========================================Menghitung cross-variogram empiris
var <- variogramST(g, data = stfdf.train,
assumeRegular = TRUE, pseudo = 1)
#write.xlsx(var, "D:/TESIS/Hasil/excels/Semivariogram empiris.xlsx")
# Ekstrak informasi dari variogramST object
gamma_values <- var$gamma # Nilai variogram (semivariance)
distances <- var$spacelag # Jarak (spatial lag)
time_lags <- var$timelag # Time lag (temporal)
# Buat plot interaktif 3D
plot_ly(x = ~distances, y = ~time_lags, z = ~gamma_values,
intensity = ~gamma_values, type = 'mesh3d',
colors = colorRamp(c("darkblue", "blue", "red", "orange", "yellow"))) %>%
layout(scene = list(
xaxis = list(title = "Jarak spasial"),
yaxis = list(title = "Lag Waktu (hari)"),
zaxis = list(title = "Gamma")
)) %>%
colorbar(title = "Gamma")
Semivariogram empiris dihasilkan menggunakan data pelatihan yang terdiri dari 23 titik lokasi selama periode 12 bulan. Komponen spasial dari semivariogram empiris menunjukkan variasi minimal dalam nilai gamma seiring dengan bertambahnya jarak spasial, yang mengindikasikan variabilitas spasial yang terbatas dalam rentang pengamatan. Sebaliknya, nilai gamma menunjukkan peningkatan yang jelas seiring dengan bertambahnya lag waktu, mencerminkan variabilitas yang lebih besar seiring waktu. Menariknya, penurunan nilai gamma diamati pada lag waktu terbesar, yang mungkin menunjukkan efek smoothing atau korelasi temporal yang berkurang pada interval yang lebih panjang.
Semivariogram empiris kemudian dimodelkan menggunakan berbagai semivariogram teoretis. Dengan menggabungkan tiga model semivariogram marginal (spherical, exponential, dan Gaussian) untuk komponen spasial, temporal, dan gabungan (spasial), total 9 model separable, 9 model product-sum, 3 model metric, 27 model sum-metric, dan 27 model simple sum-metric dievaluasi. Gambar menggambarkan bentuk semivariogram spasiotemporal yang dihasilkan dari kombinasi terbaik di setiap kategori model teoretis. Proses pemodelan ini bertujuan untuk menentukan semivariogram teoretis mana yang paling akurat menangkap pola variabilitas spasiotemporal yang diamati dalam semivariogram empiris.
head(var, 6)
## np dist gamma id timelag spacelag avgDist
## 1 0 NA NA lag0 0 days 0.000 0.0000
## 2 24 721.5611 2.119802 lag0 0 days 449.318 721.5611
## 3 12 1319.6596 8.503237 lag0 0 days 1347.954 1319.6596
## 4 60 2161.6545 4.327761 lag0 0 days 2246.590 2161.6545
## 5 72 3074.1309 2.719611 lag0 0 days 3145.226 3074.1309
## 6 96 3980.1009 6.023574 lag0 0 days 4043.862 3980.1009
tail(var, 6)
## np dist gamma id timelag spacelag avgDist
## 187 26 8604.812 2.860425 lag11 341 days 8537.042 8604.812
## 188 28 9515.440 3.537421 lag11 341 days 9435.678 9515.440
## 189 34 10284.410 17.669637 lag11 341 days 10334.314 10284.410
## 190 26 11077.128 2.925169 lag11 341 days 11232.950 11077.128
## 191 22 12132.165 30.303786 lag11 341 days 12131.586 12132.165
## 192 30 12892.806 2.563838 lag11 341 days 13030.222 12892.806
# setting lower and upper bounds
pars.l <- c(sill.s = 13.75739, range.s = 0.1*30762.69, nugget.s = 0,
sill.t = 43.72425, range.t = 31, nugget.t = 0,
sill.st = 0.1*136.6855, range.st = sqrt((0.1*30762.69^2) + (31*0)^2), nugget.st = 0,
anis = 0)
pars.u <- c(sill.s = 70.55745, range.s = 30762.69, nugget.s = 0.5*70.55745,
sill.t = 214.7959, range.t = 341, nugget.t = 0.5*214.7959,
sill.st = 136.6855, range.st = sqrt((30762.69^2) + (341*928)^2), nugget.st = 0.5*136.6855,
anis = 928)
# separable
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# separable <- vgmST("separable", space = vgm(psill = 13.75739, model = sm, range = 0.1*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), sill = 43.72425)
# separable_Vgm <- fit.StVariogram(var, separable, fit.method = 6, method="L-BFGS-B", upper = pars.u)
# MSE = c(MSE, attr(separable_Vgm, "MSE"))
# }
# }
# separable = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# MSE = as.numeric(MSE))
# separable[separable$MSE==min(separable$MSE),]
# write.xlsx(separable, "D:/TESIS/Hasil/excels/MSE Separable.xlsx")
# product-sum
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# prodSumModel <- vgmST("productSum", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), k = 0.5)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(var, prodSumModel, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# prodsum = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# MSE = as.numeric(MSE))
# prodsum[prodsum$MSE==min(prodsum$MSE),]
# write.xlsx(prodsum, "D:/TESIS/Hasil/excels/MSE Product-Sum.xlsx")
# Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_joint = c()
# MSE = c()
# for(j in daftar_model){
# model_joint = c(model_joint, j)
# metric <- vgmST("metric", joint = vgm(psill = 136.6855, model = j, range = 0.5*30762.69, nugget = 0), stAni=1)
# metric_Vgm <- fit.StVariogram(var, metric, method="L-BFGS-B", lower = pars.l, upper = pars.u)
# MSE = c(MSE, attr(metric_Vgm, "MSE"))
# }
# metric = data.frame(Joint = model_joint,
# MSE = as.numeric(MSE))
# metric[metric$MSE==min(metric$MSE),]
# write.xlsx(metric, "D:/TESIS/Hasil/excels/MSE Metric.xlsx")
# Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# for(j in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0),
# joint = vgm(psill = 136.6855, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(var, sumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# }
# summetric = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# Joint = model_joint,
# MSE = as.numeric(MSE))
# summetric
# write.xlsx(summetric, "D:/TESIS/Hasil/excels/MSE Sum Metric.xlsx")
# Simple Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# for(j in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0),
# joint = vgm(psill = 136.6855, model = j, range = sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=500, k = 0.5)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(var, SimplesumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# }
# ssmetric = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# Joint = model_joint,
# MSE = as.numeric(MSE))
# ssmetric
# write.xlsx(ssmetric, "D:/TESIS/Hasil/excels/MSE Simple-Sum Metric.xlsx")
# separable = read.xlsx("D:/TESIS/Hasil/excels/MSE Separable.xlsx")
# prodsum = read.xlsx("D:/TESIS/Hasil/excels/MSE Product-Sum.xlsx")
# metric = read.xlsx("D:/TESIS/Hasil/excels/MSE Metric.xlsx")
# summetric = read.xlsx("D:/TESIS/Hasil/excels/MSE Sum Metric.xlsx")
# ssmetric = read.xlsx("D:/TESIS/Hasil/excels/MSE Simple-Sum Metric.xlsx")
# separable[separable$MSE==min(separable$MSE),]
# prodsum=prodsum[prodsum$MSE>0,]
# prodsum[prodsum$MSE==min(prodsum$MSE, na.rm = T),]
# metric[metric$MSE==min(metric$MSE, na.rm = T),]
# summetric=summetric[summetric$MSE>0,]
# summetric[summetric$MSE==min(summetric$MSE, na.rm = T),]
# ssmetric=ssmetric[ssmetric$MSE>0,]
# ssmetric[ssmetric$MSE==min(ssmetric$MSE, na.rm = T),]
# ==============================================================MODEL TERBAIK
# separable
separable <- vgmST("separable", space = vgm(psill = 13.75739, model = 'Exp', range = 0.5*30762.69, nugget = 0),
time = vgm(psill = 109.7589, model = 'Exp', range = 217, nugget = 0), sill = 43.72425)
separable_Vgm <- fit.StVariogram(var, separable, fit.method = 6, method="L-BFGS-B", upper = pars.u)
param_sep = extractPar(separable_Vgm)
# product-sum
prodSumModel <- vgmST("productSum", space = vgm(psill = 13.75739, model = 'Gau', range = 0.5*30762.69, nugget = 0),
time = vgm(psill = 109.7589, model = 'Gau', range = 217, nugget = 0), k = 10)
prodSumModel_Vgm <- fit.StVariogram(var, prodSumModel, method = "L-BFGS-B", fit.method = 6, lower = pars.l, upper = pars.u)
param_prodsum = extractPar(prodSumModel_Vgm)
# Metric
metric <- vgmST("metric", joint = vgm(psill = 136.6855, model = 'Gau', range = 0.5*30762.69, nugget = 0), stAni=1.5)
metric_Vgm <- fit.StVariogram(var, metric, method="L-BFGS-B", fit.method = 6, lower = pars.l, upper = pars.u)
param_metric = extractPar(metric_Vgm)
# Sum-Metric
sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = 'Gau', range = 0.5*30762.69, nugget = 0),
time = vgm(psill = 109.7589, model = 'Exp', range = 217, nugget = 0),
joint = vgm(psill = 136.6855, model = 'Sph', range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
sumMetric_Vgm <- fit.StVariogram(var, sumMetric, method="L-BFGS-B", tunit="days", fit.method = 6, lower = pars.l, upper = pars.u)
param_summetric = extractPar(sumMetric_Vgm)
# Simple Sum-Metric
SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = 'Gau', range = 0.5*30762.69, nugget = 0),
time = vgm(psill = 109.7589, model = 'Gau', range = 217, nugget = 0),
joint = vgm(psill = 136.6855, model = 'Exp', range = sqrt((31000^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=0.5)
SimplesumMetric_Vgm <- fit.StVariogram(var, SimplesumMetric, method = "L-BFGS-B", fit.method = 6, lower = pars.l, upper = pars.u)
param_ssmetric = extractPar(SimplesumMetric_Vgm)
# =========================================================PLOT SEMIVARIOGRAM
sampel_plot = plot(var, wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
sep_plot = plot(var,separable_Vgm, wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
ps_plot = plot(var,prodSumModel_Vgm,wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
metric_plot = plot(var, metric_Vgm, wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
summet_plot = plot(var,sumMetric_Vgm,wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
ssmet_plot = plot(var,SimplesumMetric_Vgm,wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
sampel_plot = as.ggplot(sampel_plot)
sep_plot <- as.ggplot(sep_plot)
ps_plot <- as.ggplot(ps_plot)
metric_plot <- as.ggplot(metric_plot)
summet_plot <- as.ggplot(summet_plot)
ssmet_plot <- as.ggplot(ssmet_plot)
# Combine the plots using wrap_plots
semivar_plot <- wrap_plots(
list(sampel_plot, sep_plot, ps_plot, metric_plot, summet_plot, ssmet_plot),
ncol = 3) +
plot_layout(guides = "auto") +
plot_annotation(title = "Perbandingan Semivariogram Empiris dengan Semivariogram Teoretis")
# Print the combined plot
print(semivar_plot)
Semivariogram teoretis dengan MSE terendah menunjukkan kemampuannya untuk menjelaskan pola variabilitas yang diamati dalam semivariogram empiris. Oleh karena itu, semivariogram teoretis dengan MSE terendah akan digunakan untuk menghasilkan matriks gamma dalam persamaan prediksi cokriging. Model teoretis dengan MSE terendah adalah model sum-metric dengan nilai MSE sebesar 292,152. Selain itu, juga terlihat bahwa semivariogram sum-metric menunjukkan pola yang paling mirip dengan semivariogram empiris. Dengan demikian, model sum-metric akan digunakan untuk menghasilkan matriks gamma \(\Gamma_{11}\) untuk perhitungan bobot cokriging selanjutnya.
Proses fitting juga dilakukan terhadap cross-variogram antara variabel \(PM_{2.5}\) dengan variabel prediktor menggunakan langkah yang sama seperti pada fitting semivariogram. Fitting dilakukan menggunakan model teoretis product, product-sum, metric, sum-metric, dan simple sum-metric. Proses ini dilakukan untuk mendapatkan estimasi nilai parameter cross-variogram dengan meminimalkan selisih antara nilai empiris dengan teoretis menggunakan metode OLS. Proses ini akan menghasilkan fungsi kontinu yang akan digunakan dalam pembangkitan nilai cross-variogram.
# =========================================PM25~NO2
g_2 <- gstat(NULL, "PM25", PM25 ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
gstat("NO2", NO2 ~ 1, data = stfdf.train, fill.cross = TRUE)
# Menghitung cross-variogram empiris
cross_var_2 <- variogramST(g_2, data = stfdf.train,
assumeRegular = TRUE, pseudo = 1)
pars.l <- c(sill.s = 13.75739, range.s = 0.1*30762.69, nugget.s = 0,
sill.t = 43.72425, range.t = 31, nugget.t = 0,
sill.st = 0.1*136.6855, range.st = sqrt((0.1*30762.69^2) + (31*0)^2), nugget.st = 0,
anis = 10)
pars.u <- c(sill.s = 70.55745, range.s = 30762.69, nugget.s = 0.5*70.55745,
sill.t = 214.7959, range.t = 341, nugget.t = 0.5*214.7959,
sill.st = 136.6855, range.st = sqrt((30762.69^2) + (341*928)^2), nugget.st = 0.5*136.6855,
anis = 928)
# # separable
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# separable <- vgmST("separable", space = vgm(psill = 13.75739, model = sm, range = 0.1*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), sill = 43.72425)
# separable_Vgm <- fit.StVariogram(cross_var_2, separable, fit.method = 6, method="L-BFGS-B", upper = pars.u)
# MSE = c(MSE, attr(separable_Vgm, "MSE"))
# }
# }
# separable = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# MSE = as.numeric(MSE))
# separable[separable$MSE==min(separable$MSE),]
# # product-sum
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# prodSumModel <- vgmST("productSum", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), k = 0.5)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_2, prodSumModel, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# prodsum = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# MSE = as.numeric(MSE))
# prodsum[prodsum$MSE==min(prodsum$MSE),]
# # Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_joint = c()
# MSE = c()
# for(j in daftar_model){
# model_joint = c(model_joint, j)
# metric <- vgmST("metric", joint = vgm(psill = 136.6855, model = j, range = 0.5*30762.69, nugget = 0), stAni=1)
# metric_Vgm <- fit.StVariogram(cross_var_2, metric, method="L-BFGS-B", lower = pars.l, upper = pars.u)
# MSE = c(MSE, attr(metric_Vgm, "MSE"))
# }
# metric = data.frame(Joint = model_joint,
# MSE = as.numeric(MSE))
# metric[metric$MSE==min(metric$MSE),]
# # Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# for(j in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0),
# joint = vgm(psill = 136.6855, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_2, sumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# }
# summetric = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# Joint = model_joint,
# MSE = as.numeric(MSE))
# summetric[summetric$MSE==min(summetric$MSE),]
# # Simple Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# for(j in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0),
# joint = vgm(psill = 136.6855, model = j, range = sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=500, k = 0.5)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_2, SimplesumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# }
# ssmetric = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# Joint = model_joint,
# MSE = as.numeric(MSE))
# ssmetric[ssmetric$MSE==min(ssmetric$MSE),]
# separable$Model = c(rep("Product",9))
# separable$Joint = c(rep(" ",9))
# prodsum$Model = c(rep("Product-Sum", 9))
# prodsum$Joint = c(rep(" ", 9))
# metric$Model = c(rep("Model", 3))
# metric$Spatial = c(rep(" ", 3))
# metric$Temporal = c(rep(" ", 3))
# summetric$Model = c(rep("Sum-Metric", 27))
# ssmetric$Model = c(rep("Simple Sum-Metric", 27))
# pm_no = dplyr::bind_rows(separable, prodsum, metric, summetric, ssmetric)
# write.xlsx(pm_no, "D:/TESIS/Hasil/Cross-Variogram PM25 dan NO2.xlsx")
# Sum-Metric
sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = 'Gau', range = 0.5*30762.69, nugget = 0),
time = vgm(psill = 109.7589, model = 'Sph', range = 217, nugget = 0),
joint = vgm(psill = 136.6855, model = 'Exp', range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=1.5)
sumMetric_Vgm <- fit.StVariogram(cross_var_2, sumMetric, method="L-BFGS-B", tunit="days", fit.method = 6, lower = pars.l, upper = pars.u)
attr(sumMetric_Vgm, "MSE")
## [1] 288.9448
cv2_summetric = extractPar(sumMetric_Vgm);cv2_summetric
## sill.s range.s nugget.s sill.t range.t nugget.t
## 1.375739e+01 3.076269e+04 0.000000e+00 4.372425e+01 2.809691e+02 2.827152e+00
## sill.st range.st nugget.st anis
## 1.366855e+01 2.195445e+04 5.637747e-01 1.000000e+01
plot(cross_var_2, list(sumMetric_Vgm), all=T, wireframe=T)
Pada fitting ini, fungsi semivariogram marginal spasial menggunakan model Gaussian, sementara fungsi marginal temporal menggunakan model spherical, dan komponen joint menggunakan fungsi eksponensial. Gambar tersebut menunjukkan perbandingan cross-variogram empiris variabel konsentrasi \(PM_{2.5}\) dan konsentrasi \(NO_2\) dengan cross-variogram teoretis model sum-metric.
# =========================================PM25~Presipitasi
g_3 <- gstat(NULL, "PM25", PM25 ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
gstat("Presipitasi", Presipitasi ~ 1, data = stfdf.train, fill.cross = TRUE)
# Menghitung cross-variogram empiris
cross_var_3 <- variogramST(g_3, data = stfdf.train,
assumeRegular = TRUE, pseudo = 1)
pars.l <- c(sill.s = 914.197, range.s = 0.1*30762.69, nugget.s = 0,
sill.t = 43.72425, range.t = 31, nugget.t = 0,
sill.st = 0.1*136.6855, range.st = sqrt((0.1*30762.69^2) + (31*0)^2), nugget.st = 0,
anis = 10)
pars.u <- c(sill.s = 70.55745, range.s = 30762.69, nugget.s = 0.5*70.55745,
sill.t = 214.7959, range.t = 341, nugget.t = 0.5*214.7959,
sill.st = 136.6855, range.st = sqrt((30762.69^2) + (341*928)^2), nugget.st = 0.5*136.6855,
anis = 928)
# # separable
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# separable <- vgmST("separable", space = vgm(psill = 100, model = sm, range = 0.1*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), sill = 50)
# separable_Vgm <- fit.StVariogram(cross_var_3, separable, fit.method = 6, method="L-BFGS-B")
# MSE = c(MSE, attr(separable_Vgm, "MSE"))
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# }
# }
# separable = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# MSE = as.numeric(MSE))
# separable[separable$MSE==min(separable$MSE),]
# # product-sum
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# prodSumModel <- vgmST("productSum", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), k = 10)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_3, prodSumModel, fit.method = 6, method = "L-BFGS-B")
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# }
# }
# }
# prodsum = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# MSE = as.numeric(MSE))
# prodsum[prodsum$MSE==min(prodsum$MSE),]
# # Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_joint = c()
# MSE = c()
# for(j in daftar_model){
# model_joint = c(model_joint, j)
# metric <- vgmST("metric", joint = vgm(psill = 45000, model = j, range = 0.5*30762.69, nugget = 0), stAni=1.5)
# metric_Vgm <- fit.StVariogram(cross_var_3, metric, method="L-BFGS-B")
# MSE = c(MSE, attr(metric_Vgm, "MSE"))
# }
# metric = data.frame(Joint = model_joint,
# MSE = as.numeric(MSE))
# metric[metric$MSE==min(metric$MSE),]
# # Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# for(j in daftar_model){
# sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0),
# joint = vgm(psill = 136.6855, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=1.5)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_3, sumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# }
# }
# }
# }
# summetric = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# Joint = model_joint,
# MSE = as.numeric(MSE))
# summetric[summetric$MSE==min(summetric$MSE),]
# # Simple Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# for(j in daftar_model){
# SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0),
# joint = vgm(psill = 136.6855, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=0.5)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_3, SimplesumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# }
# }
# }
# }
# ssmetric = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# Joint = model_joint,
# MSE = as.numeric(MSE))
# ssmetric[ssmetric$MSE==min(ssmetric$MSE),]
# separable$Model = c(rep("Product", 6))
# separable$Joint = c(rep(" ", 6))
# prodsum$Model = c(rep("Product-Sum", 9))
# prodsum$Joint = c(rep(" ", 9))
# metric$Model = c(rep("Model", 3))
# metric$Spatial = c(rep(" ", 3))
# metric$Temporal = c(rep(" ", 3))
# summetric$Model = c(rep("Sum-Metric", 27))
# ssmetric$Model = c(rep("Simple Sum-Metric", 27))
# pm_ch = dplyr::bind_rows(separable, prodsum, metric, summetric, ssmetric)
# write.xlsx(pm_ch, "D:/TESIS/Hasil/Cross-Variogram PM25 dan Curah Hujan.xlsx")
# Simple Sum-Metric
SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = 'Gau', range = 0.5*30762.69, nugget = 0),
time = vgm(psill = 109.7589, model = 'Sph', range = 217, nugget = 0),
joint = vgm(psill = 136.6855, model = 'Gau', range = sqrt((31000^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=0.5)
SimplesumMetric_Vgm <- fit.StVariogram(cross_var_3, SimplesumMetric, method = "L-BFGS-B", fit.method = 6, lower = pars.l)
attr(SimplesumMetric_Vgm, "MSE")
## [1] 83946136
cv3_ssmetric = extractPar(SimplesumMetric_Vgm);cv3_ssmetric
## sill.s range.s sill.t range.t sill.st range.st
## 914.1970 659784.4230 26804.7595 238.9169 1718.8804 49970.7851
## nugget anis
## 302.3135 9728.0167
plot(cross_var_3, list(SimplesumMetric_Vgm), all=T, wireframe=T)
Dalam fitting cross-variogram antara variabel konsentrasi \(PM_{2.5}\) dan curah hujan, model simple sum-metric menghasilkan nilai MSE paling rendah dibandingkan model lainnya. Model tersebut menggunakan fungsi Gaussian pada komponen spasial, fungsi spherical pada temporal, dan fungsi Gaussian pada komponen joint.
# =========================================PM25~Kelembaban
g_4 <- gstat(NULL, "PM25", PM25 ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
gstat("Kelembaban", Kelembaban ~ 1, data = stfdf.train, fill.cross = TRUE)
# Menghitung cross-variogram empiris
cross_var_4 <- variogramST(g_4, data = stfdf.train,
assumeRegular = TRUE, pseudo = 1)
pars.l <- c(sill.s = 13.75739, range.s = 0.1*30762.69, nugget.s = 0,
sill.t = 43.72425, range.t = 31, nugget.t = 0,
sill.st = 0.1*136.6855, range.st = sqrt((0.1*30762.69^2) + (31*0)^2), nugget.st = 0,
anis = 10)
pars.u <- c(sill.s = 70.55745, range.s = 30762.69, nugget.s = 0.5*70.55745,
sill.t = 214.7959, range.t = 341, nugget.t = 0.5*214.7959,
sill.st = 136.6855, range.st = sqrt((30762.69^2) + (341*928)^2), nugget.st = 0.5*136.6855,
anis = 928)
# # separable
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# separable <- vgmST("separable", space = vgm(psill = 13.75739, model = sm, range = 0.1*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), sill = 43.72425)
# separable_Vgm <- fit.StVariogram(cross_var_4, separable, fit.method = 6, method="L-BFGS-B", upper = pars.u)
# MSE = c(MSE, attr(separable_Vgm, "MSE"))
# }
# }
# separable = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# MSE = as.numeric(MSE))
# separable[separable$MSE==min(separable$MSE),]
# # product-sum
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# prodSumModel <- vgmST("productSum", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), k = 0.5)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_4, prodSumModel, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# prodsum = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# MSE = as.numeric(MSE))
# prodsum[prodsum$MSE==min(prodsum$MSE),]
# # Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_joint = c()
# MSE = c()
# for(j in daftar_model){
# model_joint = c(model_joint, j)
# metric <- vgmST("metric", joint = vgm(psill = 136.6855, model = j, range = 0.5*30762.69, nugget = 0), stAni=1)
# metric_Vgm <- fit.StVariogram(cross_var_4, metric, method="L-BFGS-B", lower = pars.l, upper = pars.u)
# MSE = c(MSE, attr(metric_Vgm, "MSE"))
# }
# metric = data.frame(Joint = model_joint,
# MSE = as.numeric(MSE))
# metric[metric$MSE==min(metric$MSE),]
# # Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# for(j in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0),
# joint = vgm(psill = 136.6855, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_4, sumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# }
# summetric = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# Joint = model_joint,
# MSE = as.numeric(MSE))
# summetric[summetric$MSE==min(summetric$MSE),]
# # Simple Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# for(j in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0),
# joint = vgm(psill = 136.6855, model = j, range = sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=500, k = 0.5)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_4, SimplesumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# }
# ssmetric = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# Joint = model_joint,
# MSE = as.numeric(MSE))
# ssmetric[ssmetric$MSE==min(ssmetric$MSE),]
# separable$Model = c(rep("Product",9))
# separable$Joint = c(rep(" ",9))
# prodsum$Model = c(rep("Product-Sum", 9))
# prodsum$Joint = c(rep(" ", 9))
# metric$Model = c(rep("Model", 3))
# metric$Spatial = c(rep(" ", 3))
# metric$Temporal = c(rep(" ", 3))
# summetric$Model = c(rep("Sum-Metric", 27))
# ssmetric$Model = c(rep("Simple Sum-Metric", 27))
# pm_ku = dplyr::bind_rows(separable, prodsum, metric, summetric, ssmetric)
# write.xlsx(pm_ku, "D:/TESIS/Hasil/Cross-Variogram PM25 dan Kelembaban.xlsx")
# Metric
metric <- vgmST("metric", joint = vgm(psill = 136.6855, model = 'Gau', range = 0.5*30762.69, nugget = 0), stAni=1.5)
metric_Vgm <- fit.StVariogram(cross_var_4, metric, method="L-BFGS-B", fit.method = 6, lower = pars.l, upper = pars.u)
attr(metric_Vgm, "MSE")
## [1] 78.54296
cv4_metric = extractPar(metric_Vgm);cv4_metric
## sill range nugget anis
## 24.427639 22603.546735 6.308559 214.795900
plot(cross_var_4, list(metric_Vgm), all=T, wireframe=T, xlab = "Jarak Spasial", ylab = 'Lag Waktu (hari)',
zlab = 'Gamma')
Hasil fitting cross-variogram variabel konsentrasi \(PM_{2.5}\) dengan kelembaban menunjukkan bahwa model metric dengan fungsi Gaussian memberikan MSE terendah dibandingkan model lainnya
# =========================================PM25~Angin
g_5 <- gstat(NULL, "PM25", PM25 ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
gstat("Angin", Angin ~ 1, data = stfdf.train, fill.cross = TRUE)
# Menghitung cross-variogram empiris
cross_var_5 <- variogramST(g_5, data = stfdf.train,
assumeRegular = TRUE, pseudo = 1)
pars.l <- c(sill.s = 1, range.s = 0.1*30762.69, nugget.s = 0,
sill.t = 1, range.t = 31, nugget.t = 0,
sill.st = 1, range.st = sqrt((0.1*30762.69^2) + (31*0)^2), nugget.st = 0,
anis = 10)
pars.u <- c(sill.s = 3, range.s = 30762.69, nugget.s = 0.5*70.55745,
sill.t = 3, range.t = 341, nugget.t = 0.5*214.7959,
sill.st = 3, range.st = sqrt((30762.69^2) + (341*928)^2), nugget.st = 0.5*136.6855,
anis = 928)
# # separable
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# separable <- vgmST("separable", space = vgm(psill = 1, model = sm, range = 0.1*30762.69, nugget = 0),
# time = vgm(psill = 1, model = tm, range = 217, nugget = 0), sill = 1)
# separable_Vgm <- fit.StVariogram(cross_var_5, separable, fit.method = 6, method="L-BFGS-B", upper = pars.u)
# MSE = c(MSE, attr(separable_Vgm, "MSE"))
# }
# }
# separable = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# MSE = as.numeric(MSE))
# separable[separable$MSE==min(separable$MSE),]
# # product-sum
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# prodSumModel <- vgmST("productSum", space = vgm(psill = 1, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 1, model = tm, range = 217, nugget = 0), k = 0.5)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_5, prodSumModel, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# prodsum = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# MSE = as.numeric(MSE))
# prodsum[prodsum$MSE==min(prodsum$MSE),]
# # Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_joint = c()
# MSE = c()
# for(j in daftar_model){
# model_joint = c(model_joint, j)
# metric <- vgmST("metric", joint = vgm(psill = 3, model = j, range = 0.5*30762.69, nugget = 0), stAni=1)
# metric_Vgm <- fit.StVariogram(cross_var_5, metric, method="L-BFGS-B", lower = pars.l, upper = pars.u)
# MSE = c(MSE, attr(metric_Vgm, "MSE"))
# }
# metric = data.frame(Joint = model_joint,
# MSE = as.numeric(MSE))
# metric[metric$MSE==min(metric$MSE),]
# # Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# for(j in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# sumMetric <- vgmST("sumMetric", space = vgm(psill = 1, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 1, model = tm, range = 217, nugget = 0),
# joint = vgm(psill = 3, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_5, sumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# }
# summetric = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# Joint = model_joint,
# MSE = as.numeric(MSE))
# summetric[summetric$MSE==min(summetric$MSE),]
# # Simple Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
# for(tm in daftar_model){
# for(j in daftar_model){
# model_space = c(model_space, sm)
# model_temporal = c(model_temporal, tm)
# model_joint = c(model_joint, j)
# SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 1, model = sm, range = 0.5*30762.69, nugget = 0),
# time = vgm(psill = 1, model = tm, range = 217, nugget = 0),
# joint = vgm(psill = 3, model = j, range = sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=500, k = 0.5)
# # Menggunakan tryCatch untuk menangani error
# fit_result <- tryCatch({
# fit.StVariogram(cross_var_5, SimplesumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
# }, error = function(e) {
# return(NULL) # Jika ada error, kembalikan NULL
# })
#
# # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
# if (is.null(fit_result)) {
# MSE = c(MSE, 0)
# } else {
# MSE = c(MSE, attr(fit_result, "MSE"))
# }
# }
# }
# }
# ssmetric = data.frame(Spatial = model_space,
# Temporal = model_temporal,
# Joint = model_joint,
# MSE = as.numeric(MSE))
# ssmetric[ssmetric$MSE==min(ssmetric$MSE),]
# separable$Model = c(rep("Product",9))
# separable$Joint = c(rep(" ",9))
# prodsum$Model = c(rep("Product-Sum", 9))
# prodsum$Joint = c(rep(" ", 9))
# metric$Model = c(rep("Model", 3))
# metric$Spatial = c(rep(" ", 3))
# metric$Temporal = c(rep(" ", 3))
# summetric$Model = c(rep("Sum-Metric", 27))
# ssmetric$Model = c(rep("Simple Sum-Metric", 27))
# pm_angin = dplyr::bind_rows(separable, prodsum, metric, summetric, ssmetric)
# write.xlsx(pm_angin, "D:/TESIS/Hasil/Cross-Variogram PM25 dan Angin.xlsx")
# Sum-Metric
sumMetric <- vgmST("sumMetric", space = vgm(psill = 1, model = "Gau", range = 0.5*30762.69, nugget = 0),
time = vgm(psill = 1, model = "Gau", range = 217, nugget = 0),
joint = vgm(psill = 3, model = "Gau", range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
sumMetric_Vgm <- fit.StVariogram(cross_var_5, sumMetric, method="L-BFGS-B", tunit="days", fit.method = 6, lower = pars.l, upper = pars.u)
attr(sumMetric_Vgm, "MSE")
## [1] 0.1726448
cv5_summetric = extractPar(sumMetric_Vgm);cv5_summetric
## sill.s range.s nugget.s sill.t range.t nugget.t
## 1.000000e+00 3.076269e+04 0.000000e+00 1.883009e+00 3.410000e+02 2.155019e-01
## sill.st range.st nugget.st anis
## 1.000000e+00 2.684337e+04 0.000000e+00 6.152708e+01
plot(cross_var_5, list(sumMetric_Vgm), all=T, wireframe=T, xlab = "Jarak Spasial", ylab = 'Lag Waktu (hari)',
zlab = 'Gamma')
Fitting cross-variogram variabel konsentrasi PM2.5 dengan kecepatan angin menghasilkan model sum-metric sebagai model dengan MSE terendah dibandingkan dengan model lainnya. Model tersebut menggunakan fungsi Gaussian pada komponen spasial, temporal, dan joint. Gambar tersebut menunjukkan perbandingan cross-variogram empiris variabel konsentrasi PM2.5 dan kecepatan angin dengan cross-variogram teoretis model sum-metric.
sumMetric = vgmST("sumMetric", space = vgm(psill = param_summetric[1], "Gau", range = param_summetric[2], nugget = param_summetric[3]),
time = vgm(psill = param_summetric[4], "Mat", range = param_summetric[5], nugget = param_summetric[6]),
joint = vgm(psill = param_summetric[7], "Sph", range = param_summetric[8], nugget = param_summetric[9]), stAni = param_summetric[10],
temporalUnit = "days")
test.krige.ST = krigeST(PM25 ~ Presipitasi + Kelembaban + Angin + NO2,
data = stfdf.train,
newdata = stfdf.test,
modelList = sumMetric)
MAPE_test = mean(abs((test.krige.ST@data$var1.pred-stfdf.test@data$PM25)/stfdf.test@data$PM25))*100
RMSE_test = sqrt(mean((test.krige.ST@data$var1.pred-stfdf.test@data$PM25)^2))
prediksi_vs_aktual = data.frame(Prediksi = test.krige.ST@data$var1.pred, Aktual = stfdf.test@data$PM25)
e = prediksi_vs_aktual$Prediksi-prediksi_vs_aktual$Aktual
#write.xlsx(prediksi_vs_aktual, "D:/Tesis/Hasil/excels/Prediksi vs Aktual Data Test.xlsx")
cat("MAPE Test: ", MAPE_test)
## MAPE Test: 0.6588764
cat("RMSE Test: ", RMSE_test)
## RMSE Test: 0.3843981
Evaluasi dilakukan menggunakan dua ukuran kekeliruan, yaitu nilai RMSE dan MAPE. Model variogram teoretis yang dipilih memberikan prediksi yang akurat pada data uji, dengan nilai MAPE sebesar 0,659% dan RMSE sebesar 0,384.
Selanjutnya, estimasi dilakukan menggunakan model variogram dengan parameter yang telah diestimasi. Estimasi dilakukan untuk 646 grid berukuran 1 km × 1 km yang terletak secara beraturan di wilayah DKI Jakarta periode bulan Januari hingga Desember 2023.
idn.gadm = "D:/TESIS/DATA/gadm41_IDN.gpkg"
tesis.sf = st_read(idn.gadm, layer = "ADM_ADM_2")
## Reading layer `ADM_ADM_2' from data source `D:\TESIS\DATA\gadm41_IDN.gpkg' using driver `GPKG'
## Simple feature collection with 502 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.00971 ymin: -11.00761 xmax: 141.0194 ymax: 6.076941
## Geodetic CRS: WGS 84
tesis.sf = st_transform(tesis.sf, utm48s)
wilayah = tesis.sf[tesis.sf$NAME_1 == "Jakarta Raya", ]["NAME_2"]
tesis.poly = as(wilayah[-6,], "Spatial")
tesis.poly@bbox[2, 1] <- 9290000 # Ganti `new_ymin` dengan nilai y min baru
grd = SpatialPixels(SpatialPoints(makegrid(tesis.poly, cellsize = 1000)),
proj4string = proj4string(tesis.poly))
n = 12 # jumlah grid waktu/periode waktu
library(xts)
tgrd = seq(min(index(stfdf.train)), max(index(stfdf.train)), length = n)
# MENGISI GRID PREDIKSI DENGAN VERIABEL PREDIKTOR MENGGUNAKAN COKRIGING SPASIAL
# UNTUK SETIAP TITIK WAKTU
bulan = c("2023-01-01", "2023-02-01", "2023-03-01", "2023-04-01", "2023-05-01",
"2023-06-01", "2023-07-01", "2023-08-01", "2023-09-01", "2023-10-01",
"2023-11-01", "2023-12-01")
prediktor_grd = c()
for(i in bulan){
data <- stfdf.train[, i]
gstat(id = "ch", formula = Presipitasi~1, data = data) |>
gstat("rh", Kelembaban~1, data) |>
gstat("no2", NO2~1, data) |>
gstat("ws", Angin~1, data) |>
gstat(model=vgm(psill = var(data$Presipitasi), model = "Sph", range = 30000), fill.all=T) -> g
v <- variogram(g)
v.fit = fit.lmc(v, g)
plot(v, model = v.fit)
pred <- predict(v.fit, newdata = grd)
grd_pred = cbind(Presipitasi = pred@data$ch.pred,
Kelembaban = pred@data$rh.pred,
Angin = pred@data$ws.pred,
NO2 = pred@data$no2.pred)
prediktor_grd = rbind(prediktor_grd, grd_pred)
}
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
head(prediktor_grd)
## Presipitasi Kelembaban Angin NO2
## [1,] 176.0560 80.88047 4.741267 10.343881
## [2,] 176.6280 80.96468 4.734693 10.239950
## [3,] 177.2510 81.05278 4.727209 10.135030
## [4,] 177.8543 81.13615 4.719810 10.036452
## [5,] 178.4162 81.21243 4.712812 9.946155
## [6,] 178.9286 81.28097 4.706370 9.864932
# STFDF untuk new_data to predict
pred.grd = STFDF(sp = grd, time = tgrd,
data = as.data.frame(prediktor_grd))
# ==============================================================PREDIKSI FULL
tesis.krige.ST = krigeST(PM25 ~ Presipitasi + Kelembaban + Angin + NO2,
data = stfdf.train, newdata = pred.grd,
modelList = sumMetric)
colnames(tesis.krige.ST@data) = 'PM25'
min_scale = min(tesis.krige.ST@data$PM25)
max_scale = max(tesis.krige.ST@data$PM25)
# MENYIMPAN HASIL
#hasil_semua = as.data.frame(tesis.krige.ST)
#write.xlsx(hasil_semua, "D:/TESIS/Hasil/excels/Prediksi semua grid.xlsx")
# =============================================================PLOT SATU SATU
# Konversi tesis.krige.ST ke RasterLayer
raster_data <- raster(tesis.krige.ST[,2])
# Lakukan masking agar raster hanya di dalam batas wilayah
masked_raster <- mask(raster_data, tesis.poly)
# Konversi raster ke data frame
raster_df <- as.data.frame(rasterToPoints(masked_raster), xy = TRUE)
colnames(raster_df) <- c("x", "y", "PM2.5") # Sesuaikan nama kolom
# Konversi batas wilayah (tesis.poly) ke objek sf
tesis_poly_sf <- st_as_sf(tesis.poly)
tesis_poly_sf[,"NAME_ID"] = c('Jakarta Barat', 'Jakarta Pusat', 'Jakarta Selatan',
'Jakarta Timur', 'Jakarta Utara')
tesis_poly_sf[,"NAME_EN"] = c('West Jakarta', 'Central Jakarta', 'South Jakarta',
'East Jakarta', 'North Jakarta')
# ==================================================== PLOT GABUNGAN 12 BULAN
# Vektor nama bulan dari Januari hingga Desember
#bulan <- c("January", "February", "March", "April", "May", "June",
# "July", "August", "September", "October", "November", "December")
bulan = c("Januari", "Februari", "Maret", "April", "Mei", "Juni", "Juli",
"Agustus", "September", "Oktober", "November", "Desember")
# Gabungkan semua data raster ke dalam satu data frame
all_raster_df <- do.call(rbind, lapply(1:12, function(i) {
# Konversi tesis.krige.ST ke RasterLayer pada indeks tertentu
raster_data <- raster(tesis.krige.ST[, i])
# Lakukan masking agar raster hanya di dalam batas wilayah
masked_raster <- mask(raster_data, tesis.poly)
# Konversi raster ke data frame
raster_df <- as.data.frame(rasterToPoints(masked_raster), xy = TRUE)
colnames(raster_df) <- c("x", "y", "PM2.5")
# Tambahkan kolom bulan
raster_df$bulan <- bulan[i]
raster_df$time <- i
return(raster_df)
}))
# Ubah kolom 'bulan' menjadi faktor dengan urutan level yang benar
all_raster_df$bulan <- factor(all_raster_df$bulan,
levels = bulan)
# Plot dengan ggplot
peta = ggplot() +
geom_tile(data = all_raster_df, aes(x = x, y = y, fill = PM2.5)) +
geom_contour(data = all_raster_df, aes(x = x, y = y, z = PM2.5, colour = stat(level)), size = 0.2, colour = "red", alpha = 0.5) +
scale_fill_gradientn(colours = colorRampPalette(c("green4", "green", "yellow", "orangered", "red3"))(200)) +
facet_wrap(~ bulan, ncol = 4) +
theme_bw() +
ylab("") +
xlab("") +
geom_sf(data = tesis_poly_sf, fill = NA, size=0.01, colour = "black") +
geom_sf_text(data = tesis_poly_sf, aes(geometry = st_centroid(geometry), label = NAME_ID),
color = "black", size = 1.7, fontface = "bold", vjust = 1.5) + # Tambahkan label
guides(fill = guide_colorbar(title.position = "top", title.vjust = 1,
frame.colour = "black",
barwidth = 1,
barheight = 10)) +
theme(aspect.ratio = 1, legend.position = "right",
text = element_text(),
legend.direction = "vertical",
panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "grey95", linetype = "solid"),
panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
axis.text.x = element_blank(), #remove x axis labels
axis.ticks.x = element_blank(), #remove x axis ticks
axis.text.y = element_blank(), #remove y axis labels
axis.ticks.y = element_blank(), #remove y axis ticks
legend.title = element_text(size = 7),
legend.text = element_text(size = 7),
strip.text = element_text(size = 7)
) +
labs(fill = "PM2.5")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
peta
## Warning: `stat(level)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(level)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Distribusi secara spasial menunjukkan bahwa konsentrasi \(PM_{2.5}\) cenderung lebih rendah pada awal tahun (Januari hingga April) dan akhir tahun (November-Desember). Konsentrasi \(PM_{2.5}\) meningkat secara signifikan pada pertengahan tahun, yaitu selama periode Mei hingga Oktober. Pola ini berkaitan dengan faktor meteorologis, seperti meningkatnya curah hujan selama musim hujan (di akhir dan awal tahun) yang melarutkan dan menghilangkan polutan, serta curah hujan yang lebih rendah selama musim kemarau, yang membatasi proses pemurnian udara. Selain itu, kecepatan angin yang lebih rendah selama musim kemarau berkontribusi pada akumulasi polutan di area tertentu.
Hasil estimasi menunjukkan bahwa konsentrasi \(PM_{2.5}\) pada tahun 2023 berada pada rentang 15,41 µg/m3 hingga 65,88 µg/m3. Indonesia menetapkan nilai ambang batas konsentrasi PM2.5 maksimal 50 µg/m3. Wilayah-wilayah dengan konsentrasi yang melebihi nilai ambang batas sehat ditunjukkan dengan warna jingga hingga merah. Terlihat bahwa konsentrasi \(PM_{2.5}\) melebihi ambang batas pada periode-periode musim kemarau.
Konsentrasi \(PM_{2.5}\) yang tinggi umumnya teramati di Jakarta Utara, Jakarta Pusat, dan sebagian wilayah Jakarta Selatan. Di Jakarta Utara dan Jakarta Pusat, peningkatan konsentrasi \(PM_{2.5}\) terjadi karena kepadatan penduduk tinggi, aktivitas industri yang tinggi, dan emisi gas dari kendaraan bermotor. Faktor-faktor tersebut juga meningkatkan konsentrasi NO2 yang berkontribusi terhadap peningkatan \(PM_{2.5}\).
Sementara itu, tingginya konsentrasi \(PM_{2.5}\) di Jakarta Selatan dipengaruhi oleh sumber lokal dan regional. Wilayah ini berbatasan langsung dengan Bogor dan Depok. Kedua wilayah tersebut memiliki tingkat urbanisasi yang tinggi dan aktivitas transportasi yang intensif. Hal ini memungkinkan adanya polutan yang bergerak ke wilayah Jakarta Selatan. Fenomena ini terutama terlihat selama musim kemarau, yaitu ketika kondisi atmosfer (kelembaban, kecepatan angin, dan curah hujan) mendukung adanya penyebaran polutan dari wilayah sekitar.